﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace SpatialReference
{
    using System;
    public class Coordinate
    {
        public double[] coordinate;
        public Coordinate() { }
        public Coordinate(double x, double y)
        {
            coordinate = new double[] { x, y };
        }
    }
    public class TMProjected : Coordinate
    {
        public TMProjected(double x, double y)
        {
            coordinate = new double[] { x, y };
        }
        public GCoordinate Unproject(double scaleFactor, double falseEasting, double falseNorthing, double meridian, double equatorialRadius, double polarRadius)
        {
            var a = equatorialRadius;
            var b = polarRadius;
            var X = coordinate[0];
            var Y = coordinate[1];
            var m0 = scaleFactor;
            var n = (a - b) / (a + b);
            var es2 = (Math.Pow(a, 2) - Math.Pow(b, 2)) / (Math.Pow(b, 2));
            var B = (1 + n) / a / (1 + n * n / 4 + n * n * n * n / 64.0) * Y / m0;
            var b1 = 3.0 / 2.0 * n * (1 - 9.0 / 16.0 * n * n) * Math.Sin(2 * B);
            var b2 = n * n / 16.0 * (21.0 - 55.0 / 2 * n * n) * Math.Sin(4 * B);
            var b3 = 151.0 / 96.0 * n * n * n * Math.Sin(6 * B);
            var b4 = 1097.0 / 512.0 * n * n * n * n * Math.Sin(8 * B);
            var Bf = B + b1 + b2 + b3 + b4;
            var ID = Math.Floor(X / 1000000);
            var Vf = Math.Sqrt(1 + es2 * Math.Cos(Bf) * Math.Cos(Bf)); var etaf = Math.Sqrt(es2 * Math.Cos(Bf) * Math.Cos(Bf));
            var y = X - falseEasting - ID * 1000000;
            var ys = y * b / m0 / (Math.Pow(a, 2));
            var l = Math.Atan(Vf / Math.Cos(Bf) * ((Math.Exp(ys) - Math.Exp(-ys)) / 2) * (1 - etaf * etaf * etaf * etaf * ys * ys / 6 - es2 * ys * ys * ys * ys / 10));
            var lng = meridian * Math.PI / 180 + l;
            var lat = Math.Atan(Math.Tan(Bf) * Math.Cos(Vf * l) * (1 - etaf * etaf / 6 * l * l * l * l));
            lat = lat / Math.PI * 180;
            lng = lng / Math.PI * 180;
            return new GCoordinate(lat, lng, a, b);
        }
        public WGS84 UnprojectToWGS84(double scaleFactor, double falseEasting, double falseNorthing, double meridian)
        {
            var a = 6378137;
            var b = 6356752.31424518;
            var X = coordinate[0];
            var Y = coordinate[1];
            var m0 = scaleFactor;
            var n = (a - b) / (a + b);
            var es2 = (Math.Pow(a, 2) - Math.Pow(b, 2)) / (Math.Pow(b, 2));
            var B = (1 + n) / a / (1 + n * n / 4 + n * n * n * n / 64.0) * Y / m0;
            var b1 = 3.0 / 2.0 * n * (1 - 9.0 / 16.0 * n * n) * Math.Sin(2 * B);
            var b2 = n * n / 16.0 * (21.0 - 55.0 / 2 * n * n) * Math.Sin(4 * B);
            var b3 = 151.0 / 96.0 * n * n * n * Math.Sin(6 * B);
            var b4 = 1097.0 / 512.0 * n * n * n * n * Math.Sin(8 * B);
            var Bf = B + b1 + b2 + b3 + b4;
            var ID = Math.Floor(X / 1000000);
            var Vf = Math.Sqrt(1 + es2 * Math.Cos(Bf) * Math.Cos(Bf)); var etaf = Math.Sqrt(es2 * Math.Cos(Bf) * Math.Cos(Bf));
            var y = X - falseEasting - ID * 1000000;
            var ys = y * b / m0 / (Math.Pow(a, 2));
            var l = Math.Atan(Vf / Math.Cos(Bf) * ((Math.Exp(ys) - Math.Exp(-ys)) / 2) * (1 - etaf * etaf * etaf * etaf * ys * ys / 6 - es2 * ys * ys * ys * ys / 10));
            var lng = meridian * Math.PI / 180 + l;
            var lat = Math.Atan(Math.Tan(Bf) * Math.Cos(Vf * l) * (1 - etaf * etaf / 6 * l * l * l * l));
            lat = lat / Math.PI * 180;
            lng = lng / Math.PI * 180;
            return new WGS84(lat, lng);
        }
    }
    public class GCoordinate : Coordinate
    {
        public double a;
        public double b;
        public double e2;
        public GCoordinate() { }
        public GCoordinate(double[] coordinate, double equatorialRadius, double inverseFlattening)
        {
            this.coordinate =  coordinate;
            this.a = equatorialRadius;
            this.b = equatorialRadius * (1 - 1 / inverseFlattening);
        }
        public GCoordinate(double latitude, double longitude, double equatorialRadius, double polarRadius)
        {
            coordinate = new double[] { latitude, longitude };
            this.a = equatorialRadius;
            this.b = polarRadius;
        }
        public GCoordinate HelmertTransform(double px, double py, double pz, double ks, double rw, double rf, double rk)
        {
            var lat = coordinate[0] * Math.PI / 180;
            var lng = coordinate[1] * Math.PI / 180;
            e2 = (Math.Pow(a, 2) - Math.Pow(b, 2)) / Math.Pow(a, 2);
            var v = a / Math.Sqrt(1 - e2 * (Math.Sin(lat) * Math.Sin(lat)));
            var X = v * Math.Cos(lat) * Math.Cos(lng);
            var Y = v * Math.Cos(lat) * Math.Sin(lng);
            var Z = (v * (1 - e2)) * Math.Sin(lat);
            var c1 = Math.Cos(rw);
            var c2 = Math.Cos(rf);
            var c3 = Math.Cos(rk);
            var s1 = Math.Sin(rw);
            var s2 = Math.Sin(rf);
            var s3 = Math.Sin(rk);
            var D11 = c2 * c3;
            var D21 = -c2 * s3;
            var D31 = Math.Sin(rf);
            var D12 = s1 * s2 * c3 + c1 * s3;
            var D22 = -s1 * s2 * s3 + c1 * c3;
            var D32 = -s1 * c2;
            var D13 = -c1 * s2 * c3 + s1 * s3;
            var D23 = c1 * s2 * s3 + s1 * c3;
            var D33 = c1 * c2;
            var X1 = px + ks * (D11 * X + D12 * Y + D13 * Z);
            var Y1 = py + ks * (D21 * X + D22 * Y + D23 * Z);
            var Z1 = pz + ks * (D31 * X + D32 * Y + D33 * Z);
            lng = Math.Atan2(Y1, X1);
            var lat0 = Math.Atan2(Z1, Math.Sqrt(X1 * X1 + Y1 * Y1));
            e2 = (Math.Pow(a, 2) - Math.Pow(b, 2)) / Math.Pow(a, 2);
            while (Math.Abs(lat - lat0) > 0.0000000001)
            {
                var No = a / Math.Sqrt(1 - e2 * Math.Sin(lat0) * Math.Sin(lat0));
                var h = Math.Sqrt(X1 * X1 + Y1 * Y1) / Math.Cos(lat0) - No;
                lat = lat0;
                lat0 = Math.Atan(Z1 / Math.Sqrt(X1 * X1 + Y1 * Y1) * (1 / (1 - e2 * No / (No + h))));
            }
            return new GCoordinate(lat / Math.PI * 180, lng / Math.PI * 180, a, b);
        }
        public TMProjected TransverseMercartor(double scaleFactor, double falseEasting, double falseNorthing, double meridian)
        {
            var lat = coordinate[0] * Math.PI / 180;
            var lng = coordinate[1] * Math.PI / 180;
            lng = lng - meridian * Math.PI / 180;
            var m0 = scaleFactor;
            var es2 = (Math.Pow(a, 2) - Math.Pow(b, 2)) / (Math.Pow(b, 2));
            var V = Math.Sqrt(1 + es2 * Math.Cos(lat) * Math.Cos(lat));
            var eta = Math.Sqrt(es2 * Math.Cos(lat) * Math.Cos(lat));
            var Bf = Math.Atan(Math.Tan(lat) / Math.Cos(V * lng) * (1 + eta * eta / 6 * (1 - 3 * Math.Sin(lat) * Math.Sin(lat)) * lng * lng * lng * lng));
            var Vf = Math.Sqrt(1 + es2 * Math.Cos(Bf) * Math.Cos(Bf));
            var etaf = Math.Sqrt(es2 * Math.Cos(Bf) * Math.Cos(Bf));
            var n = (a - b) / (a + b);
            var r1 = (1 + n * n / 4 + n * n * n * n / 64) * Bf;
            var r2 = 3.0 / 2.0 * n * (1 - n * n / 8) * Math.Sin(2 * Bf);
            var r3 = 15.0 / 16.0 * n * n * (1 - n * n / 4) * Math.Sin(4 * Bf);
            var r4 = 35.0 / 48.0 * n * n * n * Math.Sin(6 * Bf);
            var r5 = 315.0 / 512.0 * n * n * n * n * Math.Sin(8 * Bf);
            var Northing = a / (1 + n) * (r1 - r2 + r3 - r4 + r5) * m0 - 0.001 + falseNorthing;
            var ys = Math.Tan(lng) * Math.Cos(Bf) / Vf * (1 + etaf * etaf * lng * lng * Math.Cos(Bf) * Math.Cos(Bf) * (etaf * etaf / 6 + lng * lng / 10));
            ys = Math.Log(ys + Math.Sqrt(ys * ys + 1));
            var Easting = m0 * Math.Pow(a, 2) / b * ys + falseEasting;
            return new TMProjected(Easting, Northing);
        }
    }
    public class WGS84 : GCoordinate
    {
        public WGS84(double latitude, double longitude)
        {
            coordinate = new double[] { latitude, longitude };
            a = 6378137;
            b = 6356752.31424518;
        }
    }
    public class GreekGrid
    {
        public static double[] WorldToGreek(double[] world)
        {
            return new WGS84(world[0], world[1]).HelmertTransform(199.723, -74.03, -246.018, 1, 0, 0, 0).TransverseMercartor(0.9996, 500000, 0, 24).coordinate;
        }
        public static double[] GreekToWorld(double[] greek)
        {
            var coord = new TMProjected(greek[0], greek[1]).UnprojectToWGS84(0.9996, 500000, 0, 24).HelmertTransform(-199.723, 74.03, 246.018, 1, 0, 0, 0).coordinate;
            return new double[] { coord[0], coord[1] };
        }
    }
}



